﻿module KMean
open MathNet.Numerics
open MathNet.Numerics.FSharp
open MathNet.Numerics.LinearAlgebra.Double
open System
open System.Collections.Generic
open Types
open Utils

let r_min (xs : List<DenseVectorWrapper>) (n : int) (mu : DenseVectorWrapper list) =
  let xn = xs.[n]
  
  mu  |> Seq.map    (fun p -> 
                      let dist = xn.DenseVector - p.DenseVector
                      let sqdist = dist * dist
                      sqdist, p)
      |> Seq.minBy  (fun (dist, p) -> dist)
      |> (fun (_,p) -> p.Index)


let mutable clusters_changed = false  
let rec populate_clusters (mu_k_list : DenseVectorWrapper list) (mu_k_points_list : HashSet<DenseVectorWrapper> list) (input : List<DenseVectorWrapper>) =
  for observed_point in input do
    let observerd_mu_index = r_min input observed_point.Index  mu_k_list
    if observerd_mu_index <> observed_point.PrevCluster.Value then clusters_changed <- true
    observed_point.PrevCluster := observerd_mu_index
    mu_k_points_list.[observerd_mu_index].Add(observed_point) |> ignore  

let avg_mu (mu_k : DenseVectorWrapper) (mu_k_points : HashSet<DenseVectorWrapper>) =
  let mu_k_points_length = mu_k_points.Count
  let densevector_acc = DenseVector(mu_k.DenseVector.Count)
  if mu_k_points_length = 0 then densevector_acc
  else
    let mu_k_points = mu_k_points |> Seq.toList
    mu_k_points |> List.fold(fun s t -> t.DenseVector + s) densevector_acc    
                |> (fun tot -> tot / float mu_k_points_length)
  

let optimization_of_mu (mu_k_list : DenseVectorWrapper list) (mu_k_points_list : HashSet<DenseVectorWrapper> list) (input : List<DenseVectorWrapper>) =  
  for mu_k_points_list in mu_k_points_list do mu_k_points_list.Clear()
  clusters_changed <- false
  do populate_clusters mu_k_list mu_k_points_list input
  let mu_k_list = mu_k_list 
                        |> List.rev
                        |> List.fold(fun s t -> DenseVectorWrapper((avg_mu t mu_k_points_list.[t.Index]), t.Index, -1, ref t.Index) :: s) []
  clusters_changed |> not, mu_k_list

let make_random_vector (input : List<DenseVectorWrapper>) =
  input.[seed.Next(input.Count)].DenseVector

let k_mean (k : int) (input : List<DenseVectorWrapper>)=

  let mu_k_list =
    [for k in [0..k-1] do   
      //mu with supervised builder
//      let k_point = input |> Seq.filter(fun p -> p.RealNumber = k) |> Seq.toList
//      let hash_set = HashSet(k_point)
//      let vector = avg_mu k_point.[0] hash_set
//      yield DenseVectorWrapper(vector , k, -1, k |> ref)]
      
//      //mu without supervised builder
      yield DenseVectorWrapper((make_random_vector input), k, -1, k |> ref)]

//      //enrico compare
//      yield DenseVectorWrapper(input.[k * 100].DenseVector, k, -1, k |> ref)]

  let mu_k_points_list = [for _ in mu_k_list do yield HashSet<DenseVectorWrapper>()]

  //do populate_clusters mu_k_list mu_k_points_list input

  let rec find_mu i mu_k_list (prev_mu_k_list : DenseVectorWrapper list) =
    
    let res,new_mu_k_list = optimization_of_mu mu_k_list mu_k_points_list input

    do graphics.Clear(white_color)

    for mu_k, prev_mu_k in List.zip mu_k_list prev_mu_k_list do
      let prev_mu_k_matrix = prev_mu_k.toMatrix16x16
      let mu_k_matrix = mu_k.toMatrix16x16
      do draw_point mu_k_matrix prev_mu_k_matrix mu_k.Index

    f.Update()
    Console.WriteLine(string i)
    if res then mu_k_list, i
    else find_mu (i + 1) new_mu_k_list mu_k_list 

  let res = find_mu 0 mu_k_list (mu_k_list.Tail @ [mu_k_list.Head])

  for mu_k, prev_mu_k in List.zip (fst res) (fst res) do
    let prev_mu_k_matrix = prev_mu_k.toMatrix16x16
    let mu_k_matrix = mu_k.toMatrix16x16
    do draw_point mu_k_matrix prev_mu_k_matrix mu_k.Index

  //compute R
  let mutable a = 0.
  let mutable b = 0.
  for p1 in input do
    for p2 in input do
      if p1.Index <> p2.Index then 
        if !p1.PrevCluster = !p2.PrevCluster && p1.RealNumber = p2.RealNumber then
          a <- a + 1.
        elif !p1.PrevCluster <> !p2.PrevCluster && p1.RealNumber <> p2.RealNumber then
          b <- b + 1.
  let r = (2. * (a + b)) / (float input.Count * (float input.Count - 1.))
  Console.WriteLine("R = " + string r)
  res